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ABSTRACT 

Size differences of « 20% between red (metal-rich) and blue (metal-poor) sub- 
populations of globular clusters have been observed, generating an ongoing debate 
as to weather these originate from projection effects or the difference in metallicity. 
We present direct TV-body simulations of metal-rich and metal-poor stellar popula- 
tions evolved to study the effects of metallicity on cluster evolution. The models start 
with N = 100 000 stars and include primordial binaries. We also take metallicity 
dependent stellar evolution and an external tidal field into account. We find no signifi- 
cant difference for the half-mass radii of those models, indicating that the clusters are 
structurally similar. However, utilizing observational tools to fit hali-light (or effec- 
tive) radii confirms that metallicity effects related to stellar evolution combined with 
dynamical effects such as mass segregation produce an apparent size difference of 17% 
on average. The metallicity effect on the overall cluster luminosity also leads to higher 
mass-to-light ratios for metal-rich clusters. 

Key words: globular clusters: general - galaxies: star clusters: general - stars: mass- 
loss - stars: luminosity function, mass function - methods: iV-body simulations 



1 INTRODUCTION 

Globular clusters (GCs) are substantial components of 
galaxies and found in populations of up to thousands in giant 
elliptical galaxies l|Peng et al.ll201lh . The Milky- Way (MW) 
hosts a GC population of 157 confirmed clusters (|Harrid 
1 19961 . 20 10 edition), with n ew clusters still being discov- 
ered (e.g. iMinniti et al]|201ll 'l. These clusters live within the 
bulge as well as the halo of the Galaxy and can - in contrast 
to star clusters beyond the Local Group - easily be resolved 
in ground-based observations. In general, the GC systems 
of galaxies tend to ap pear in two sub-populations: a blue 
and a red component (|Zinnl ll985). Although the metallic- 
ity cannot be inferred directly fro m the cluster co lour due 
to the age-metallicity degeneracy (|Worthevl [l99l h it has 
been well accepted that the blue clusters are metal-poor, 
whereas the r ed ones are metal-rich. Bo th sub-populations 
are old (e.g. iMarfn-Franch et al.l 120091 '). with a trend for 
the red clusters to be more centrally concentrated within 
their host galaxy's potential than the ir blue counterparts 
jKinmanll 19591 : iBrodie fc Straderl r2006 , also see Fig. [lj. 

The ability of the Hubble Space Telescope to par- 
tially resolve globular clusters even beyond the Local Group 
has lead to the finding that i) GCs have mean half-light 



E-mail: asippel@astro.swin.edu.au 



radii r-^i = 3 pc (J ordan et al.l 120051 ) and ii) red clusters 
are on avera ge ~ 17 — 30% smaller than their metal poor 
counterparts jKundu fc Whitmorelll998l ; Ijordan et a"l]|2005l : 
IWoodlev fc G6mej|2010r T Several explanations for this phe- 
nomena have been prop osed: projection effe cts and the in- 
fluence of the tidal field (|Larsen et al .1 120011 or a combined 
effect of mass segregation and the dependence of main- 
seque nce lifetimes on metallicity |jordanll2004l ; I Jordan et all 
2005). Whether either of those effects are dominating or a 
combination of both can only be investigated through direct 
star cluster simulations where three-dimensional galactocen- 
tric distances are known and stellar evolution is included in 
the dynamical evolution of the cluster. 

The effects of metallicity on the evolution of a single 
star manifests itself as a different rate of stellar evolution, 
which is accompanied by a different mass-loss rate and hence 
ultimately affects the stars lifetime and remnant mass (see 
Section f2|. In general, low metallicity stars evolve faster 
along th e main sequence th an their high metallicity coun- 
terparts <|Hurlev et aLl koOO). For a bound system such as 
a star cluster, the increased mass-loss rate can lead to a 
lower cluster mass and hence a lower escape velocity and. 
This in turn can produce a stronger increase in radius for 
the metal-poor cluster. At later stages, this might also lead 
to postponed core-collapse for the low metallicity cluster. 
Both effects could lead to a larger measured cluster size. A 



2 Anna C. Sippel, Jarrod R. Hurley, Juan P. Madrid, William E. Harris 



preliminary study a long these lines has been carried out by 
iHurlev et al] \2QQ4 ) for open clusters. They showed that an 
increased escape rate for the metal-rich clusters owing to 
earlier core-collapse acts to cancel these effects resulting in 
only a 10% difference in cluster lifetime for metal-poor ver- 
sus metal-rich cases - within the statistical noise of fluctuat- 
ing results from one simulation to another. However, several 
aspects of our new simulations differ from this preliminary 
study. Among those are an adjusted binary fraction for GCs 
and an improved tidal field. Most importantly we also use a 
higher initial number of stars Ni, bringing the iV-body mod- 
els into the GC regime. This ultimately leads to an increase 
in cluster lifetime and hence not necessarily core-collapse or 
depletion of stars within a Hubble time. 

In this work, we make use of a set of star cluste r simula- 
tions evolv ed with the direct iV-body code NB0DY6 l)Aarsethl 
1 19991 . [2003) to study the effects of metallicity on star cluster 
dynamics, evolution and size (i.e. effective radius) to an- 
swer the question if metallicity alone could reproduce the 
observed size difference. We measure the sizes of these clus- 
ters along their evolutionary track with methods used both 
in observation s and theory. 

Recently IDownind HmJ) has published a set of Monte- 
Carlo models exploring the origin of the observed size dif- 
ference between metal-rich and metal-poor GCs, which pro- 
vides an excellent comparis on for our work. This f ollows on 
from the iV-body models of ISchulman et all l|2012T h who in- 
vestigated the evolution of half-m ass radius with m etallicity 
in small-iV clusters. Similarly to IDownind (|2012h . we shall 
be careful to make a distinction between the actual size of 
a star cluster, represented by the half-mass radius (which 
we shall denote as r 50 %, i.e. the 50% Lagrangian radius), 
and the observationally determined size (the half-light or 
effective radius, r c ff). 

This paper is structured as follows. We introduce the 
differences in stellar evolution depending on metallicity in 
the next Section. In Section [3] we describe our simulation 
method and the models we have chosen to evolve. In Section 
[4] we analyze the evolution of cluster mass, binary fraction, 
luminosity, half-light radius and mass-to-light ratio which is 
followed by discussion and conclusions. 



2 METALLICITY EFFECTS ON STELLAR 
AND STAR CLUSTER EVOLUTION 

The main sequence (MS) lifetime of a single star depends 
mainly on its mass (and hence luminosity), but al so on its 
chemi cal composition: the metallicity Z (or [Fe /H] ) . Iciavtonl 
( 1968) shows that the MS lifetime can be represented as: 
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where m and L are a star's mass and luminosity and X the 
hydrogen fraction. A star's mass at given luminosity scales 
as: 



(2) 



where ko is the central opacity and [t, the mean molecular 
weight. The hydrogen fraction X and helium abundance Y 
can be set as a function of metallicity according to: 
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Figure 1. Size and gal actocentric di stance of the MW GC pop- 
ulation (compiled from Harris 2010). Blue circles are used for 
metal-poor and red squares for metal-rich clusters, the distinction 
is made at [Fe/H]= —1.1. The solid black like denotes the size- 
distance relation ry ss *J d sc from Ivan den Bergh et al ] 4l99ll y 
Metal-poor clusters tend to have larger galactocentric distances 
as well as larger sizes (half-light radii). The models used for this 
study are evolved at a galactocentric distance of 8.5 kpc, marked 
by the vertical dotted line. 



X = 0.76 - 3Z 



(3) 



Y = 0.24 + 2Z (4) 

as m iPols et all (|1998T ). If Z is decreased, X increases while 

Y decreases slightly, leading to a marginally lower mean 
molecular weight: 

o 

(5) 



1 + 3X + 0.5Y 



To first order, it can be assu med that the central opacity 
is proportional to Z: ko oc Z |Clavtonlll96sl ). Using this, in 
combination with Eqs. [1] and [2] we find: 
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with kJJ 2 oc Z ' 2 being the dominant term in this equation. 
A lower opacity implies less resistance for escaping photons 
from the hydrogen burning core and hence a higher luminos- 
ity and therefore a shorter lifetim e (see also Tab le [p. For 
an extended discussion we refer to Iciavtonl (|l968T l. 

In the iV-body models, we e volve stars ac c ordin g to the 
stellar evolution prescriptions of lHurlev et al. (|2000l), which 
are based on the detailed models of lPols et alj l| 19981 ). These 
prescriptions are accurate for a wide range of metallicities 
and cover all phases of stellar evolution. This means stars are 
evolved from the zero-age main sequence up to and includ- 
ing the remnant phases: white dwarfs (WDs), neutron stars 
(NSs) and black holes (BHs). If necessary, the stellar evolu- 
tionary track evolves via the giant branch, core helium burn- 
ing and t hermally pul s ating asymptotic giant branch. As 
shown bv lHurlev et al.l |2000i ). the difference in MS lifetime 
is most prominent for low-mass stars and steadily decreases 
towards higher mass stars until M m 8 Mq , where the high 
metallicity stars begin to have a shorter MS lifetime, al- 
though only marginally (and noting that model uncertain- 
ties are more prevalent at higher masses). This implies, that 
for clusters of the same age, the mass of the most massive MS 
star (and hence MS turnoff mass m-ro) is higher in a high-Z 
cluster. Examples for m-ro are given in Table [1] It is not 
only the MS lifetime that is altered by the metallicity, but 
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Table 1. Main sequence lifetimes for stars with different metallicities. Metallicity Z and [Fe/H] are in the first two columns, followed by 
the hydrogen (X) and helium (Y) mass fraction (Eq. \3\ and |4j . Even though the mean molecular weight fi (Eq. [5j in column 5 is barely 
affected by the metallicity, different relative MS lifetimes *ms (column 6) are caused by a change in o pacity for different metallicities 
according to Eq. [6]. The expected MS lifetimes up to the Hertzsprung Gap according to lHurlev et al,l bOOOT ) for stars with initially 3, 
i.5 and 0.8 Mq are given in the next three columns, followed by the MS turnover mass jjito in columns iO and if at ages of if and 
i2Gyr. We note that stars with Z = O.OOi and Z = O.OOO f evolve in a s i milar fashion compared to the metal rich case - hence Z = 0.001 
is also a metal-poor case. This has already been noted bv lHurlev et al. I J2004h . as well as the fact that stars and clusters with Z = 0.01 
evolve similar to solar metallicity Z = 0.02. 
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0.0001 


-2.3 


0.76 


0.24 


0.59 


0.25 


0.29 Gyr 


1.6 Gyr 


13.5 Gyr 


0.84 Mq 


0.83 Mq 


0.001 


-1.3 


0.76 


0.25 


0.59 


0.40 


0.29 Gyr 


1.7 Gyr 


14 Gyr 


0.85 Mq 


0.83 M e 


0.01 


-0.3 


0.74 


0.26 


0.60 


0.60 


0.35 Gyr 


2.4 Gyr 


21.7 Gyr 


0.95 Mq 


0.91 Mq 



also the remnant mass. For initial masses less than 50 Mq 
our models give a maximum black hole mass tubu ~ 28 Mq 
for metal -poor stars versus msH ~ 12 Mq for metal rich pro- 
genitors l|Belczvnski et al1l2006t ). This trend is the same for 
all remnants: a 2 Mq progenitor with Z — 0.0001 will end 
life as a WD of mass m = 0.84 Mq, while a metal-rich coun- 
terpart with Z — 0.01 will have a WD mass of m = 0.66 Mq. 
This occurs after w 0.9 and 1.4 Gyr, respectively. Hence the 
remnant mass in a metal-poor cluster is always expected to 
be higher (see also Table 0. 

Since there is no strong evidence for an ex- 
plicit meta llicity dependen ce of the mass-los s rate 

of giants (llben fe Renzinil 1 19831 : ICarraro et al. i urn 
ISchroder fc Cunt j I2005T T generally mass-loss from the 
envelope during the giant branch phase and beyond 
is implemented according to Reimer's law (formula of 
iKudritzki fc Reimers|[l978T ): 

LR _1 in\ 

m oc Mq yr . (7) 

An implicit metallicity dependence exists as the evolution 
of the radius R and L depend on the mean molecular weight 
and hence Z, as mentioned earlier (e.g. Eq. HJ. Exceptions 
apply for very massive stars, e.g. luminous blue variables 
with luminosity L > 4000 Lq , where the mass- loss is mod- 
eled according to: 
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This is Eq. 2 from [Nieuwen huii zen fc de Jagerl (I1990T) but 
modified by the factor Z - 5 (|Kudritzki et al.l I1989T ). Note 
that mass-loss can also occur as a result of mass transfer 
- having ultimately the same effect of moving a star along 
the MS towards lower effective temperature and hence lower 
luminosity. 

2.1 Stellar evolution of an entire population 

To quantify the effects of stellar evolution on a non- 
dynamical population, we evolv e 105 000 stars to gether 
through stellar evolution alone |Hurlev et al.l [2000). This 
means that dynamical effects such as the influence of the 
galactic tidal field as well as the intrinsic ./V-body evolu- 
tion within the cluster are ignored. The set-up of the initial 
masses of this population is identical to our iV-body mod- 
els introduced in Section O where the dynamical evolution 
is fully incorporated. In Fig. [2] the mass, luminosity and 
mass-to-light ratio evolution of this model is illustrated for 



the three metallicities Z = 0.01, Z = 0.001 and Z = 0.0001. 
At the Hubble time, « 30% of the initial stellar mass is lost 
purely due to stellar evolution and only « 50% of the ini- 
tial mass in MS stars is stil l remaining (in agreement with 
iBaumgardt fc Makind l2003h ■ The overall mass of the low- 
Z population stays higher throughout, while the mass con- 
tained in MS stars is always higher in the high-Z population, 
as expected due to the higher MS turnoff mass. The luminos- 
ity (actually calculated as the V-band luminosity) drops by 
an order of magnitude within the first « 2 Gyr and roughly 
another magnitude over the next 10 Gyr of evolution. We see 
that even though a high-i? cluster will have a higher 7tito, 
the luminosity of a metal-poor cluster remains 1.5 — 2 times 
higher throughout the entire evolution - based on stellar 
evolution alone. This implies that the increased brightness 
of \ow-Z stars is outweighing the higher number of MS stars 
in the high-Z case. As expected from the evolution of mass 
and luminosity in this non-dynamical model, the mass-to- 
light ratio M/L is predicted to be higher by nearly a factor 
two for a metal-rich cluster. In a dynamically evolved model 
with a tidal field, the mass-to light ratios are likely to be 
modified as preferential ly low-mass stars are lost fr om the 
outskirts of the cluster (jBaumgardt fc Makino 2 003] ) . Low- 
mass MS stars are faint and have a high mass-to-light ratio. 
We will compare Fig. [5] to our dynamical models in Section 
1431 



2.2 



Size: projection effects vs. internal dynamics 

For the MW. Ivan den Bergh et all l|l99ll ) found that the GC 
half-light radius ru can be related to the galactocentric dis- 
r d^ (see Fig. [TJ . As the MW is the only 



tance d s 
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galaxy where three-dimensional galactocentric distances are 
available, one has to rely on projected distances for extra- 
galactic GC systems. Studies of extragalactic GC systems 
have shown that red and blue clusters are found to have 
different spatial distributions within the potential of their 
host galaxy: red clusters are distributed closer to the cen- 
tre of the galaxy and subject to a stronger influence of th e 
tidal field than the blue clusters (Brodic & Stradcr 2006). 
A siz e difference ranging f rom 17% (| Jordan et al.l 120051 ) to 
30% l|Woodlev et alJIioOcf ) for the blue and red population 
has been found in numerous studies. Several scenarios have 
been proposed for the origin of this effect: projection effects 
and/or the effect of stellar evolution in combination with 
mass segregation, which we describe below. In addition, the 
possibility of different initi al conditions during cluster for- 
mation have been proposed (|Harrisll200g| l as well as different 
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Figure 2. Population of 105 000 stars, corresponding to the set-up for our ./V-body models, evolved with only stellar evolution (i.e. no 
dynamical interactions). The metallicities are Z = 0.0001 (dashed blue), Z = 0.001 (solid green) and Z = 0.01 (dashed-dotted red). Left 
panel: The upper set of lines is the mass (scaled by the initial mass) of the entire population of stars (including all remnants) while the 
lower set of lines only includes stars on the main sequence. This illustrates that the mass contained in MS stars is always higher for 
metal-rich clusters and that the overall mass evolution depends critically on the treatment of remnants. Middle panel: Corresponding 
luminosity evolution (in units of Lq). In this case the treatment of remnants is not crucial. Note that even though the metal-rich model 
contains more stars on the MS, the metal-poor models have a higher luminosity. This is in agreement with Eq. [2] and implies that the 
higher luminosity of low-Z stars is outweighing the fact that the metal-poor clusters have a lower MS turnoff mass at any given time 
than the metal-rich cluster (see Table [TJ . Note that both metal-poor cases are expected to evolve in a similar fashion (also see Table [T{ , 
however some variation is caused depending on the number of bright stars at any given time. Right panel: Resulting mass to light ratio 
M/L (in units of Mq/Lq). As expected from the luminosity evolution, the mass- to-light ratio is higher for the metal-rich model. 



initial mass functions for metal-poor or metal-rich clusters 
(jStrader et alJl20QgL IpqTlh . 

lLarsen &i Brodiel (|2003T ) found that projection effects 
may account for the observed size difference of red and blue 
GCs, if the GC distribution flattens out near the centre of 
the galaxy (e.g. King profile) and there is a steep relation 
between cluster size r and galactocentric distance d gc . How- 
ever, this is not the case for either more centrally peaked dis - 



tributions or shallower r — d sc relations. Spitler et alj (|2006l ) 



find in agreement with lLarsen &: Brodk ( 20031 ) that proiec- 
tion effects could explain the observed size difference in the 
Sombrero galaxy. A size gradient for GCs is found for small 
but not large (projected ) galactocentri c distances. 

In contrast to this, IJordanl \2004 l has found that the 
combined effects of mass segregation and MS lifetime lead 
to a size difference of low vs. high metallicity clusters. Un- 
der the assumption that the average half-mass radius does 
not depend on metallicity, the observed light-profiles were 
modeled with Michie-King multi-mass models and stellar 
isochrones leading to the result that a size difference of 
the observed magnitude arises naturally, with the metal- 
rich model having a half-light radius ~ 14% smaller than its 
metal-poor counterpart. The reasoning for this originates 
from the different speed in stellar evolution of stars with 
different Z implying that the light profile of a high Z clus- 
ter can appear more concentrated. Unfortunately, in this ap- 
proach the interplay between stel lar dynamics a nd evolution 
was not considered. We note that I Jordan! l)2004l ) assumes the 
average half- mass radius to be independent of [Fe/H] - an 
assumption pointing to a universality in the formation and 

evolution of GCs. 

As part of the ACS Virgo Cluster Survey l|C6te et al.l 

|2004 ) the sizes of thousands of globular clusters belong- 
ing to 100 early ty pe ellipticals in Virgo were measured 
I Jordan et al.ll2005l ). They find in agreement with previous 
studies that the average half-light radius depends on the 
color of the GCs, with red GCs being « 17% smaller than 
their blue counterparts. This size difference was proposed to 



originate from the effects of mass-segregatio n and metallic - 
ity, hence intrinsic cluster mechanisms as in I Jordan! (|2004l ). 

The arguments given above show that it is necessary to 
know the three dimensional galactocentric distance of GCs 
to their host galaxy to fully understand and disentangle the 
influence of the environment and metallicity on GC evo- 
lution. To be able to distinguish between those effects, we 
focus on metal-poor and metal-rich clusters at the same loca- 
tion, i.e. where both coexist. In the MW, 16 GCs are located 
in the region between 7 < d gc < 10 kpc with a mean size of 
rhi = 3.95 pc. Four of those are metal-rich ([Fe/H]> 1.1) and 
12 metal-poor. Thus we chose a galactocentric distance of 
dgc = 8.5 kpc for our models. 



SIMULATION METHOD & CHOICE OF 
PARAMETERS 



We use the direct N-body code NB0DY6 dAarsethll 19991 , 12003T ) 
to construct and evolve our models. This state-of-the art 
A^-body code takes advantage of the possibility to carry 
out such simulations on a graphics processing unit (GPU) 
coupled together with conv entional central processing units 
jNitadori fc Aarsethl I2012T I. The simulations were carried 
out on Tesla S1070 graphics cards at Swinburne University. 

We use a K roupa initial mass function (IMF: 
iKroupa et ail 1 19931 ) within the limits 0.1 to 50A/q to pop- 
ulate our cluster model with stars. The beginning t = 
for the simulation corresponds to the zero-age MS and no 
gas is included in the models. The simulations start with 
Ni — 100 000 stellar systems, including a primordial binary 
frequency of 5% (see Section [3. If) . These stars are initially 
distributed following a Plummer density profile 



p(r) = 



3M 



1 + 



{rsc) 



1 -5/2 



(9) 



i|Plummeilll91ll ; lAarseth et~ailll974l ) where M is the cluster 
mass and R 3C is a scale radius (see below). As the Plummer 
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profile formally extends out to infinite radius, a cut-off at 
ten times the half-mass radi us is applied to avoid rare cases 
of stars at large distances (|Aarsethi r2003). The individual 
initial positions and velocities are then assigned such that 
the cluster is in virial equilibrium. 

The cluster is subject to a constant, MW-like tidal 
field consisting of three co mponents: a point-mass bulge, 
an extended smooth disc ijMivamoto fc Nagail [l975l. and 
a dark matter halo. We use M b = 1.5 x 1O 1O M and 
M d = 5 x 10 10 Mp for bulge and disc mass, respectively 
l|Xue et al.l I2008T ). The scale parameters for the Miyamoto 
disc are a — 4kpc (disc scale length) and 6 = 0.5 kpc 
(galactic thickness). Formally the disk extends to infinity 
but with this choice of parameters the strength has dropped 
to less than 0.1% of the central value at a distance of 
40 kpc. The dark ma tter halo follow s a logarithmic profile 
$ oc Vq ln(d 2 + 6 2 ) ' 5 (|Aarsethll20"0^ 'l. Here d is the distance 
from the galactic centre at any given time, and b is con- 
strained such that the combined mass of the bulge, disk and 
halo give an orbital velocity of vo = 220km/s at a galacto- 
centric distance of d gc = 8.5 kpc. 

As mentioned earlier, we choose to place our clusters in 
an orbit at d gc ~ 8.5 kpc to match an environment where 
red and blue clusters coexist within the MW (see Fig. [T]). 
The orbit is inclined w 22 deg to the galactic disc reaching a 
maximum height of z » 3 kpc above the galactic plane. The 
apogalacticon is 8.8 and perigalacticon 8.2 kpc with orbital 
period of » 0.2 Gyr (see Fig.[3|. We chose a mid eccentricity 
to not start with extreme cases. The inclination results in a 
maximum z = 3 kpc, w hich is typical for many MW clusters 
l|Dauphole et al . 1996). During the lifetime of a cluster, stars 
are naturally lost due to dynamical relaxation, evolution and 
disc-shocking events. The tidal radius of a cluster in the 
Milky Way potential described above can be approximated 
as: 

n ^{2w) (10) 

l|Kiipper et al.ll2010h . where f2 is the angular velocity of the 
cluster orbit and G is the gravitational constant. Calculated 
at apogalacticon gives a n = 52 pc, which we take as our ini- 
tial value. This is adjusted as the cluster evolves according 
to the factor M 1//3 . Stars are only removed from the clus- 
ter once their di stance from the clu ster centre exceeds twice 
the tidal radius. iGieles et al.l (|201lh have expressed the im- 
pact of the galactic tidal field on a cluster by quantifying a 
boundary M\i m < 1O 5 M0 x 4kpc/R gc below which clusters 
are tidally affected, whilst more massive clusters are tidally 
unaffected. The clusters in this study fall below this limit 
and hence are tidally limited. 

Within the framework of NB0DY6, the only remaining 
parameter is the scale radius i? sc , which sets the initial clus- 
ter size or density and acts as a conversion factor between 
physical and A^-body units. It is an ongoing debate as to how 
extended GCs are when they are born. Recently, it has been 
pointed out that GCs could be the remnants of much bigger 
stel lar structures such as the nuclei of accreted dwarf ga lax- 
ies l|Freemanlll993l : lBoker|[2008l : iForbes fc Bridges! l201Ch . In 
general, shortly after stellar nuclear fusion is ignited within 
a proto-cluster, the cluster it is expected to increase it's size 
as the remaining gas not incorporated into stars during star 
formation is ejected from the cluster. So far, globular clus- 
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Figure 3. Mass in the outskirts of the cluster (black solid line) 
and height z above the galactic plane (dashed grey line). M ou t is 
defined as the mass between one and two tidal radii. The galactic 
disc corresponds to z = 0. Equivalent behaviour is observed for all 
models and metallicities. Approximately 30 Mq are lost at every 
disc crossing. 

ter sizes at this early stage cannot be determined through 
observations. We choose R sc = 8, corresponding to an ini- 
tial three-dimensional half-mass radius of r 5 o% ~ 6.2 pc. The 
half-mass radius evolves to w 7 pc at the Hubble time, but is 
a three dimensional quantity and hence a smaller half-mass 
radius by 25% would be expected when me asuring projected 
radii in two dimensions (|Fleck et al.ll2006h . This places our 
models within the size range of observed clusters in the MW 
at d gc ~ 8.5 kpc (se e Fig. [!} as well as in the large and small 
Magellanic Clouds l|Mackev et al.ll2008h . 

3.1 Binary fraction 

All models used in this study are evolved with the same 
number of initial stellar systems, Ni = 100 000, incorporat- 
ing a binary fraction 

bf = AT N \ r =0.05. (11) 

This translates into N s = 95 000 single stars and Nt = 5000 
binary systems, therefore 105 000 stars in total. Some of 
these primordial binary systems may be disrupted early on, 
while new binaries form during the cluster evolution due to 
two- or three-body interactions. Open cluster studies found 
in the lit erature are usually evolved with binary frac tions of 
0.2-0.5 l|Hurlev et al.ll2004l , l2005l ; lTrenti et al.ll2007h . as ob- 
servations find higher binary fractions in th ese clusters (e.g. 
iMonteomerv et al.|[l993l : iRicher et al]|l998l for M67). Much 
lower binary fractions are observed in GCs: iMilone et al.l 
|2012l) have measured the binary fractions of 59 GCs in 
the MW and commonly find values around bf « 0.05. Bi- 
nary systems in GC models have proven to be important 
from a dynamical point of view: even a small binary frac- 
tion in the core can be sufficient to heat the cl uster core 
enough to postpone co re-collapse significantly l|Hut et al.l 
ll992l : lHeggie et alJlioOrj 1 ). 

Within NB0DY6, standard binary evolution is t reated 
according to the binary algorithm of iHurlev et al.l (|2002f ) 
where circularization of eccentric orbits as well as angular 
momentum loss mechanisms are modeled. Wind accretion 
from one binary component to the other is possible as well 
as mass transfer when either star fills its Roche lobe. Sta- 
ble hierarchi cal three- and four-body s ystems are detected 
and evolved l|Mardling fc Aar scth 2001), with single-binary 
and binary-binary encounters followed directly. This allows 
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Figure 4. Mass- loss rate for a cluster model dynamically evolved (including stellar evolution) - compared to Fig. [2] without dynamical 
interactions. Left Panel: The solid line denotes the overall cluster mass (note that in contrast to Fig. [2]the high-Z population is no longer 
less massive than the low-Z populations). The dashed line is the mass contained in MS stars and the dotted line is the mass contribution 
from WDs. Middle panel: Corresponding luminosity evolution. Although most apparent up to ?s 7 Gyr, the metal-poor clusters stay 
more luminous throughout the entire evolution. As exected the overall luminosity is lower than in the non-dynamical case (Fig. Right 
panel: Resulting mass-to-light ratio. Within the noise the values are equivalent to Fig. [2] Hence the overall dynamical evolution has an 
impact on the mass of the cluster, but little effect on M/L. 



for the replacement of one member of a binary by an in- 
coming star, formati on of binaries in few-body encounters 
and direct collisions (|Kochaneklll993 ). often leading to the 
formation of exotic stars such as blue stragglers. Nearby 
stars can perturb binary s ystems and cause chaotic orbits 
jMardling fc Aarseth|[200ll ). 



3.2 Treatment of remnants 

Neutron stars are assumed to be subject to a velocity kick 
arising from asymmetries during their formation through 
core-collapse supernovae, with observations of NSs indicat- 
ing a vast range of velocities from several up to hundreds of 
km/s. Such velocities are easily in excess of a typical GC es- 
cape velocity and, in combination with observations of sub- 
stantial NS populations in GCs, is known as the neutron- 
star retentio n problem JPfahl et al.l |2002| ). Indeed, X-ray 
sources (e.g. IWoodlev et al.ll2008l in the case of NGC 5128) 
and milli-second pulsars ( Bogdanov et al.ll201ll in the core 
of NGC 6626) indicate that NSs and BHs are common and 
even BH-BH binaries may exist. In A-body simulations, sev- 
eral different methods to assign velo city kicks to NS or BH 
remna nts have been used in the past. lBaumgardt fc Makind 
(2003) simply retain all NSs. With their IMF not reaching 
masses higher than 15 Mp, the number of NSs is not exces- 
sive and no BHs form. lMackev et al.l (|2007| ) retain all stellar- 
mass remnant BHs whilst using an IMF up to 100 Mq . In 
contrast to this, Zonoozi et alj (|201 lh retain no NSs or BHs. 
iHurlev fc Mackevl (|2010T )~ use a Gaussian velocity kick distri- 
bution peaked at « 190 km/s for both NSs and BHs, where 
the formation of a BH-BH binary is later observed to post- 
pone core-collapse. 

In this study, we adopt an intermediate approach by 
choosing V], at random from a fiat kick distribution in 
the range — 100 km/s and assigning this to NSs and 
BHs at their birth. Because of the low escape velocity 
v e = y/2GM/r w 4.7 km/s at the half-mass radius, or 
v e ~ 2.8 km/s at the tidal radius (both at a cluster age 
of 500 Myr), this reproduces a retention fraction of ~ 5% 



20021 ). We use the same al gorithm to assign kic k 
velocities to BHs at their formation l|Repetto et al.|]2012l ). 



We note that the metallicity influences the mass of the rem- 
nants. In our model, the maximum BH mass is ~ 30 Mq for 
metal-poor p rogenitors and ~ 1 Mq for their metal- rich 
counterparts (|Hurlev et al.ll2000l ; iBelczvnski et al.ll2010| ). 



3.3 MODELS 

We evolve three sets of models a), b) and c) with identical 
set-up apart from the random number seed for the initial 
particle distributions. Each set consists of three models with 
metallicities Z = 0.0001, Z = 0.001 and Z = 0.01 (see 
Table [2}, i.e. low, intermediate and high metallicity. GCs 
in the MW are found w ithin the metallicity range —2.37 < 
[Fe/H] < l|Harrisl ll996). The intermediate metallicity case 
Z = 0.001 of this study already corresponds to a metal- 
poor cluster in the MW (and also other galaxies). The low- 
Z case Z = 0.0001 is an example from the metal-poor end 
of the metallicity distribution. We expect these two low- 
metallicity clusters to exhibit similar evolution to each other 
(e.g. the MS turnoff masses agree fairly well: see Table[TJ but 
distinct from t he high- meta ll icity case. This has previously 
been noted bv lHurlev et al.l |2004). All models are evolved 
up to 14 Gyr, whil e we concentrate ou r analysis at typical 
GC age of 12 Gyr l|Hansen et al.ll2007l ). 



4 EVOLUTION 

During cluster evolution, stars are lost in three ways: i) an 
increase of velocity during two-body encounters (evapora- 
tion) or ejection following three or four-body encounters, ii) 
velocity kicks owing to SN explosions and iii) tidal strip- 
ping and disc shocking (i.e. the influence of the tidal field of 
the host galaxy). These three effects cannot be completely 
disentangled: the former two might bring a star close to 
or even beyond the tidal boundary r t (Eq. 1 10 p such that 
when crossing the galactic disc, stars are easily removed 
from the system. The periodicity of this event is « 100 Myr 
and causes the number of stars within the cluster envelope 
between one and two tidal radii to continually fluctuate be- 
tween 180 — 240 Mq , with ~ 30 Mq lost each time the disc is 
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Table 2. Mctallicitics and initial masses for all models at t = and various parameters at 12 Gyr: mass M, number of stars N, binary 
fraction bf (Eg. lilt , number of MS stars Nms as well as the total mass contained in MS stars Mms> mass locked in WDs Mwd and the 
number of NSs and BHs, all at 12 Gyr. The variation in the initial cluster mass arises from the difference in random seed when drawing 
stars from the IMF. Note the consistently higher WD mass for metal-poor clusters: WD masses are higher for same-mass progenitors 
and in addition more stars have already turned off the MS into WDs. 
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[Fe/H] 
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M (12 Gyr) 


N 


bf 




M MS 


M WD 


N NS 


JVbh 


a) 


0.0001 


-2.3 


6.43 x 10 4 


1.57 x 10 4 


35699 


0.0597 


27626 


9.48 x 10 3 


6.14 x 10 3 


18 


2 




0.001 


-1.3 


6.43 x 10 4 


1.51 x 10 4 


34787 


0.0595 


27025 


9.40 x 10 3 


5.51 x 10 3 


31 


2 




0.01 


-0.3 


6.43 x 10 4 


1.45 x 10 4 


35680 


0.0604 


27975 


1.03 x 10 4 


4.44 x 10 3 


21 


2 


b) 


0.0001 


-2.3 


6.42 x 10 4 


1.56 x 10 4 


35958 


0.0606 


27948 


9.57 x 10 3 


6.10 x 10 3 


25 


2 




0.001 


-1.3 


6.42 x 10 4 


1.50 x 10 4 


34654 


0.0595 


27036 


9.39 x 10 3 


5.42 x 10 3 


23 
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0.01 


-0.3 


6.42 x 10 4 


1.51 x 10 4 


35017 


0.0615 


28229 


1.04 x 10 4 


4.38 x 10 3 


16 
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c) 


0.0001 


-2.3 


6.36 x 10 4 


1.59 x 10 4 


36149 


0.0585 


28065 


9.66 x 10 3 


6.13 x 10 3 


26 







0.001 


-1.3 


6.36 x 10 4 


1.52 x 10 4 


34804 


0.0606 


27015 


9.43 x 10 3 


5.50 x 10 3 


24 


2 




0.01 


-0.3 


6.36 x 10 4 


1.49 x 10 4 


34596 


0.0628 


27896 


1.03 x 10 4 


4.32 x 10 3 


27 


1 



3000 - 




Figure 5. Mass function of the dynamically evolved stellar pop- 
ulation at 12 Gyr for Z = 0.01 (left panel), Z = 0.001 (middle 
panel) and Z = 0.0001 (right panel). Here we focus on model set 
b) but the behaviour is similar for all sets. The grey area is the 
entire population of stars, the thin black line the remaining stars 
on the main sequence and the dashed line the contribution of 
white dwarfs (peaked at rs 0.6 Mq). For stars with M < 0.5Mq 
the population is made entirely out of MS stars. For metal-poor 
clusters, the MS turnoff is noticeably smaller (see also Table [TJ. 
The number of NSs and BHs is insignificant compared to MS stars 
and WDs. 



crossed (see Fig. [3j) . The evolution of a star cluster is linked 
to the two-body relaxation time: 



0.14JV 



ln(0.4iV) V GM 



50% 



(12) 



JSpitzer fc Hard |l97ll : iBinnev fc Tremainei I200&1 see also 



iHurlev et all 12001 



For our models, the relaxation time is 
highest at w 2.5 Gyr when t r h corresponds approximately to 
the cluster lifetime at that point. The relaxation time then 
decreases to roughly one Gyr at 12 Gyr of cluster age. The 
half-mass relaxation timescale is not significantly affected 
by the metallicity: variations up to 10% occur. This means 
that the clusters of different metallicity are dynamically of 
similar age, which is not the case from a stellar evolution 
point of view. As can be seen in Table [2] all three metallic- 
ity models have similar mass and number of stars at 12 Gyr, 
whereas the distribution of mass among MS stars and rem- 
nants differs (the metal-poor cluster containing more mass 
in remnants). 

Fig. U is a reproduction of Fig. [5] but now for our full 
./V-body models. The three panels are again mass, luminosity 
and M/L. In Fig. [2] we only considered the mass-loss owing 



to stellar evolution, causing « 40% mass-loss up to 12 Gyr. 
For Fig.[4]the dynamical interactions are taken into account, 
resulting in an additional mass-loss of the same order, leav- 
ing a cluster mass of w 25% after 12 Gyr. Hence the three 
effects i)-iii) mentioned above are together responsible for 
approximately half the mass-loss of the cluster, while stellar 
evolution alone is responsible for the rest of the mass-loss. 
It can also be seen in Fig. [4] that nearly 40% of the mass 
at 12 Gyr is contained in WDs. All stars are split into their 
relevant stellar populations in Fig. [5] for further illustration, 
also at a cluster age of 12 Gyr. While stars below 0.5 Mq 
exclusively are on the main sequence, the contribution of 
WDs is significant for higher masses, causing a second peak 
in the mass function at » 0.6 Mq. 

There is always more (luminous) mass contained in MS 
stars in the metal-rich cluster, which is expected from the 
higher MS turnoff mass (see Table [T] and Fig. [5]), while the 
low-Z clusters are more luminous in spite of this. See Section 
14.51 for further details of the evolution of luminosity and 
mass-to-light ratio. Even though the MS turnoff is higher 
for metal-rich clusters, the number of MS stars is not always 
the highest (see Table O column 7). The fluctuation of Nms 
is mainly due to the fact that the number of low-mass stars 
with m < 0.2 Mq varies depending on metallicity: the mass 
in those stars is typically 5 — 10% higher for Z = 0.0001 than 
for Z — 0.001 or Z — 0.01. While statistical noise may be 
responsible for some of the fluctuations, the lowest- Z cluster 
also has the highest mass and hence a slightly higher escape 
velocity. 



4.1 Binary systems 

The binary fraction of initially 0.05 slightly increases to 
« 0.06 at 12 Gyr for any metallicity or model (see Table 
[2j, where the binary fraction for the high-Z model is al- 
ways slightly higher than for low metallicities. While some 
of the initial systems may easily disr upt, others form during 
few-b ody encounters. Hard binaries jHeggic 1975: lHut et al.l 
I1992D have been shown to successfully halt core-collapse over 
large periods of time and BH binary systems in particular 
can heat the core substantially. With an initial mass function 
up to 50 Mq and the inclusion of stellar evolution, black hole 
remnants will form early on in the cluster evolution. While 
most BHs get ejected almost immediately (e.g. Section [3]2]), 
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remaining BHs will sink towards the centre of the cluster ow- 
ing to mass segregation. While doing so, they may become 
part of a binary or triple system, breaking up a previously 
existing binary system. Once BHs are part of binary sys- 
tems, BH-BH binaries can easily form in a further encounter 
through exchange interactions. All BH-BH binaries in this 
stydy are dynamical binaries, having formed through such 
few-body interactions. This also means each component in 
a BH binary is the remnant of a high-mass MS progenitor 
that was either a single star or born in a binary system that 
later disrupted. 

4.2 Cluster size 

Owing to the cumulative effects of mass-loss, two-body re- 
laxation and the influence of the tidal field, the models 
are expected to go through an initial expansion, followed 
by contraction. In Fig. [6] this is shown by means of the 
three-dimensional half-mass radius r 50 %. We find no metal- 
licity dependence on the half-mass radius. Moving further 
inwards, we look at the 10% Lagrangian radius r 10 % (mid- 
dle panels of Fig. [5} and the iV-body core radius r c (bottom 
panel). Small differences in r 10 % are evident for the differ- 
ent models, noting that this inner radius is susceptible to the 
actions of highly energetic binaries in the core, even so, the 
evolution of r 10 % remains fairly steady. The iV-body core 
radius r c is similar in size to the 10% Lagrangian radius, 
however we see that r c is heavily fluctuating when BHs, 
BH-BH binaries or otherwise energetic binary systems are 
present (all of which are more likely to reside in the central 
regions owing to mass segregation). 

The iV-body core radius is not to be confused with 
an observational King-core radius, as the iV-body core ra- 
dius is a density weighted mean distance to the cluster 
centre (not taking luminosity into account). In the pro- 
cedure of calculating r a> the mean density (in terms of 
mass) of the six neighbo ring stars is calculated for each star 
l|Casertano fc Hutlll985f ). introducing a large bias towards 
stars in the neighborhood of BHs: a BH might be up to 
28 Mq , a binary BH up to twice as much, while a MS is less 
massive than two solar masses after one Gyr of cluster evo- 
lution. We note that the A^-body core radius is consistently 
less fluctuating at high metallicity than in the lower metal- 
licity cases. This is not a sampling effect. Instead, it results 
from remnant masses being lower for the high-Z population. 
With BH masses only up to 10 Mq, the density contrast 
around stars will be less steep. BH-BH binaries can mimic 
core-collapse (Fig. [SJ) when indeed just a subsystem of stars 
is responsible for this effect. 

In addition, peaks in the core radius can (but don't 
have to be) closely correlated with highly-energetic binary 
systems. As an example, the drop in r c for the \ow-Z model 
b) in Figure [6] (middle panel) at 2 Gyr is caused by a short- 
period binary composed of two carbon-oxygen white dwarfs 
of masses 0.7 and 0.8 M . At t = 1.85 Gyr, the two WDs 
merged and the product was subsequently ejected from the 
cluster. The maximum binding energy before coalescence is 
141 Mq/AU. This is followed by another dip in r c at 2.6 Gyr, 
when the core radius shrinks to 0.72 pc. At this point, more 
than 50% of the core- mass is contained in BHs and a BH-BH 
binary forms. 

In the \aw-Z model of set a), the iV-body core radius 



drops by more than factor of two to 1.4 pc at 5 Gyr. This is 
caused by a chain of reactions involving four remnant BHs 
(out of ten present at that time). The masses of the four 
BHs are 27, 26, 14 and 11 Mq, respectively. Initially, the 
least massive BH is ejected from this four-body subsystem, 
and leaves the cluster. The remaining three form a short- 
lived triple-system which ends with a BH-BH binary and 
a single BH being ejected from the core as a result from 
enhanced velocities obtained in the interaction. This implies 
that four of the most massive components are lost from the 
core within a time frame of only 40Myr. 

We conclude that the metallicity has no effect on the 
half-mass radius or other scaling parameters based on cluster 
mass. However as Figs.[2]and[4]already indicate - the metal- 
licity influences the overall luminosity of GCs with high-Z 
clusters being fainter than metal-poor clusters. To explore 
this possibility in more det ail, we meas ure the half-light of 
effective radius r e g by fitting[King ( 1966) models to our clus- 
ters - analogous to sizes are measured from observations. We 
illustrate this method in Section T4. 31 

4.3 Surface brightness and half light radii 

Among other properties, the output of NB0DY6 incorporates 
the mass, luminosity and radius for each star. This means 
effective temperatures can easily be calculated and this data 
can be cr oss convolved with stellar atmosphere model cal- 
culations (|Kuruczlll979l l to obtain Johnson V-band magni- 
tudes. 

We project this data in a two dimensional image and 
slightly smooth it with a Gaussian filter (see Fig. for an 
example of a cluster at the age of 13 Gyr). This means the 
light of each star is conserved, but not contained within one 
single pixel, which implies that the starlight can be divided 
between consecutive bins when creating a surface brightness 
profile, which is crucial in cases of very bright stars. For each 
model, at each snapshot three such images are obtained by 
using the degree of freedom to project in either the x, y or 
z direction (in theory mult iple projections are possible, see 
iNovola fc Baumgarddl^llT ) and a surface brightness profile 
is obtained separately for each projected snapshot (Fig. [7]). 
For s implicity, w e assume a background of zero. We chose 
tofit lKmg|<|l966h models as they have shown to be a robust 
solution to fit GCs. Another option would be I Wilson] (|l975l ) 
models, ha ving a greater sensitivity in the outer regions of 
the cluster (|McLaughlin et ai1l2008l ). However, in this work 
we are not investigating tidal fluctuations but the overall 
cluster evolution, which the King models are well suited for. 
Since there is no analytical solution for the surface density of 
this model, a grid of mo del fits has to be pre-calc ulated. We 
utilize the gridf it code (|McLaughlin et al"1l2008T l where this 
has been done. Each snapshot is fitted three times according 
to the three different projections along the x-, y- and z-axes, 
as illustrated in Fig. \7\ Obvious bad fits are rejected from 
further analysis (note that no bright stars have been masked 
for fitting). For each given time, the final effective radius is 
the mean along all three projections. 

The result is plotted in Fig.[8]over the entire evolution of 
the cluster. Similar to the half-mass radius, an initial expan- 
sion when mass-loss is dominated by stellar evolution winds 
from massive stars in the core is followed by a contraction 
when the mass-loss is dominated from the cluster bound- 
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Figure 6. Evolution of lagrangian radii and core radius with models a) on the left, b) in the middle panel and c) on the right. Top: 
half-mass radius. Slight size differences between the models occur, however this is not primarely related to metallicity. Size differences are 
originating from few-body encounters and high-energy binary systems in the centre of the cluster and are enhanced in the 10% lagragian 
radius (middle). Bottom: The TV-body core radius r c . Short-term effects on the core radius are often linked to high-energetic binaries or 
the presence of BHs in the core, which can severely impact the evolution of r c . 




Figure 7. Example fits for a cluster at the age of 13 Gyr. The three panels denote the same cluster at the same time, projected along 
the x—, y— and z— axis. The corresponding snapshot is printed above. Each snapshot is fitted individually. The measured data points 
for the surface brightness profile are denoted by black sqares with poisson error, the black line is the gridf it King66 fit. The resulting 
effective radius r e g is indicated by the red dotted line. 
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Figure 8. Half-light or effective radius r e g from lKfn3 l|l966t ) model fits using gridf it (McLaughlin ct al. 2008T). In the top panels, the 
overall evolution of the half-light radius is indicated for all sets of models: a) on the left, b) in the middle and c) on the right. Of greatest 
interest is the data at late times, which are highlighted below. Average cluster sizes for each metallicity are calculated for the intervalls 
10.25 — 11 Gyr, 11 — 11.75 Gyr, 11.75 — 12.5 Gyr and 12.5 — 13.25 Gyr, using blue squares for the low-Z, green diamonds for intermediate 
and red circles for the high-Z case. It is clearly seen, that the metal-poor cluster snapshots (blue) have a larger observed half-light radius 
than the metal-rich (red) snapshots. The average sizes are summarized in Table [3] 



ary. Yet there are differences in comparison to the half-mass 
radius: Firstly, r e s is approximately half as large as the half- 
mass radius. As r$Q% is a three dimensional quantity, r e ff is 
expected to be only 3/4 as large simply due to projection 
effects. A size difference further to this implies that the lu- 
minosity alters the measured cluster size. Secondly, there is 
a clear effect of the metallicity on the r e s evolution of the 
clusters: the metal-poor clusters are consistently observed 
to be larger than their metal-rich counterparts. 

Also in Fig. [8] we highlight the time window of 10 — 
13 Gyr which is of most significance for old GCs. The data 
is averaged over St — 750 Myr windows: tio = 10.25 — 11 
Gyr, tn = 11 - 11.75 Gyr, tvx = 11.75 - 12.5 Gyr and 
ti3 = 12.5 — 13.25 Gyr. The results are summarized in Ta- 
bleland combined give an overall size difference of « 17% 
between red and blue GCs. If split into sets, the difference 
is 19, 21 and 10% for sets a), b) and c), respectively. This 
result implies that the observed size difference between the 
metal-poor and metal-rich GC sub-populations can (at least 
partly) be explained by the effects of metallicity. 

4.4 Origin of the size difference and influence of 
remnants 

We observe no size difference with metallicity for the clusters 
when measuring the size by means of the mass distribution, 



Table 3. Average cluster sizes measured for all sets for the 
intervals t 10 = 10.25 - 11 Gyr, tu = 11 - 11.75 Gyr, ti2 = 
11.75 - 12.5 Gyr and t 13 = 12.5 - 13.25 Gyr. In the bottom line 
the size difference Ar = r\> — r T is given for the corresponding 
time interval, where is the average cluster size observed for 
blue, metal-poor and r r for red, metal-rich clusters. The overall 
size difference for all ages is 17%. 

tio ill *12 *13 

Z = 0.01 4.30 pc 4.08 pc 3.85 pc 3.82 pc 
Z = 0.001 4.82 pc 4.81 pc 4.61 pc 4.39 pc 
Z = 0.0001 5.01 pc 4.75 pc 4.64 pc 4.31 pc 

A r 16.5% 16.4% 20.5% 12.6% 

e.g. half-mass radius. This indicates that the clusters are 
structurally identical, and different mass-loss rates depend- 
ing on metallicity are not causing the cluster size to change 
appreciably. Also, the overall mass and mass segregation are 
not largely affected by metallicity: a higher MS turnoff mass 
for the metal-rich cluster is compensated by a lower remnant 
mass, two effects almost canceling each other out. In Fig.[5]i) 
we show the typical radial profile of the average stellar mass 
for the three different metallicities at a late age. The mod- 
els are in good agreement, showing no significant variation 
with Z. However we find size differences of up to 20% when 
measuring the cluster size by means of the stellar luminosity. 
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Table 4. Luminosity L and mass- to-light ratio M/L for stars with 
different masses and metallicities. For given mass, the luminosity 
increases with metallicity, causing M/L to decrease. 
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Figure 10. Evolution of the luminosity contained withing the 
10% Lagrangian radius normalized by the total luminosity at that 
time, for model set b). The high-Z cluster (dotted red) has a 
higher concentration of light within r 10 % than the metal-poor 
models. 



The reason for this is two- fold. Firstly, less massive remnants 
in the high-Z cluster free more space in the core for MS and 
giant stars, i.e. luminous matter, steepening the luminosity 
profile in the central regions. This is evident in Fig. [5] ii) 
which plots the radial profile of the average luminosity per 
radial region. The second factor can also be clearly seen in 
the same figure: even though low-Z clusters have a lower MS 
turnoff mass, the luminosity of MS stars of identical masses 
is higher in the \ow-Z case. This results in the low-Z clusters 
appearing brighter beyond the centre, with the differences 
beyond two parsecs being significant in relation to the error- 
bars, as shown in Fig.[5]ii). Combined, these effects result in 
a larger cluster appearance for the metal-poor clusters. To 
reinforce this we show in Fig. [10] the luminosity within the 
10% Lagrangian radius normalized by the total luminosity, 
as a function of time. Here we see that the metal-rich cluster 
consistently has a greater central concentration of luminous 
matter. 

The fact that low-Z stars are brighter for a given mass 
than their metal-rich counterparts, will be the case inde- 
pendent of a different treatment for NSs and BHs. However, 
different NS and BH abundances might affect the surface 
brightness profile by altering the central concentration of 
luminous stars. We have evolved an additional set of models 
where NSs and BHs receive a larger kick at formation, re- 
sulting in neither sub-population being present in the cluster 
after a few hundred Myr of cluster evolution (with the ex- 
ception of the rare case that a NS may form via a WD- WD 
merger). In contrast to the previous models that contain 
NSs and BHs, this causes the luminosity profiles for differ- 
ent metallicity clusters to be nearly identical (see the far 
right panel of Fig.[9|. This is no surprise: the remnant mass 
depends on metallicity and removing the remnants erases 
some of the metallicity eff ects. This i s in e xcellent agree- 
ment with the findings by iDownind <|2012h . where signif- 
icant half-light radii differences are measured with Monte 
Carl o models (utilizing the same stellar evolution prescrip- 
tion |llurie^et|al][200o|) only when BHs are retained in the 
cluster. While o ur model clusters are smaller than those of 
IDownind |20L2T ). and we only retain a few BHs compared 
to hundreds in their study, we find the same effect already 
with very few BHs present, with a contribution also from 
the NSs that are present. 



4.5 Mass-to-light ratio 

In Section \2. II we have already mentioned the mass-to- light 
ratio M/L for a stellar population evolved purely with stel- 
lar evolution, but no dynamical interaction (see Fig. [2] right 
panel). The higher overall luminosity for metal-poor pop- 
ulations implies a lower M/L ratio: the mass-to-light ra- 
tio increases with increasing metallic ity. The same trend 
has previously been observed by e.g. lAnders et al.l (|2009l ) 
w here GALEV models wer e com puted based on the models 
of iBaumgardt fe Makinol ||2003h . In Fig. g] we repeated the 
same analysis as in Fig. f2] but now for our JV-body models. 
We chose model set b) as an illustrative case, but all three 
sets are equivalent. The evolution of mass for all metallicities 
is nearly identical (Fig. [3]), whereas the metal-poor cluster 
has a slightly higher overall mass while the metal rich clus- 
ter has a slightly higher MS mass. The overall luminosity is 
evolving in a similar fashion as in the non-dynamical model, 
but a factor of two lower owing to the loss of stars. Metal- 
licity differences in L are obvious especially for t < 6 Gyr, 
but continue up to 13 Gyr. The dynamical evolution in- 
troduces selective effects on the evolution of M/L as low- 
mass main sequence st ars are preferentially lost fr om the 
outskirts of the cluster (|Baumgardt fe M akino 2003). Those 
low- mass stars have a high M/L. White dwarfs also have 
relatively low average mass compared to stars in the cen- 
tral regions. Thus they are candidates to be lost and have 
a mass-to-light ratio approaching infinity. As a general rule, 
losing a low-mass MS star or a white dwarf will decrease the 
mass-to- light ratio (see Table 0]). There is an additional ef- 
fect arising from metallicity differences to consider: for any 
given mass at a certain time, the luminosity of the metal- 
poor star will be higher than for a metal-rich star and hence 
the \aw-Z star will have a lower M/L. This implies that es- 
caping metal-rich stars will cause a larger decrease of M/L. 
In other words: the mass-to-light ratio will be more affected 
by the loss of low-mass stars in a hig h-Z cluster. While this 
is in a gree ment with the mode ls by IBaumgardt fc Makinol 
|2003r i and lAnders et al.l l|200Gt) . it is in disagreement with 
the o bserved mass-to-light rati o s of metal-rich clust ers in 
M31 (|Strader et al J 120091 lioTTI l. IStrader et al.l (|201lh have 
suggested different initial mass functions for red GCs, which 
has not been tested here. 



5 DISCUSSION AND CONCLUSIONS 

We have measured the sizes of GC models with different 
metallicity, evolved with the direct A-body code NB0DY6. 
All clusters start their evolution with 105 000 stars and a 
mass of « 6 x 10 4 Mq. We find no size differences with 
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Figure 9. Radial profiles of the average stellar mass (panel i) and average luminosity (panel ii)) for different metallicities. This is for 
model b) at an age of 11.5 Gyr, averaged over ten consecutive snapshots (covering about 130 Myr). The shaded regions indicate the 
errors involved, calculated as the standard deviation from the mean within those ten snapshots. In calculating the average mass all stars 
and remnants are taken into account, while for the average luminosity only stars not yet in the remnant phase are taken into account 
(e.g. only luminous stars). There is a general trend for the luminosity distribution to be steeper in the high metallicity case. However 
beyond the core of the cluster, the metal-poor cluster has a higher average luminosity. The panels iii) and iv) are a repeat of the panels 
on the left, but for a set of models without NS or BH remnants. While the overall evolution for these models is similar to the other 
models in this study, we do not observe a significant size difference. 



metallicity when measuring sizes by means of the half-mass 
radius or other mass-weighted radii, with the exception that 
lower remnant masses for high-Z stars cause the TV-body 
core radius to fluctuate less. This indicates, that there is 
no structural difference between clusters of low and high 
metallicity. Even though the mass-loss rates of low-Z stars 
are higher, especially in the initial stages of evolution, a con- 
sequently lower escape velocity and higher average remnant 
mass cancels this effect, leading to no overall size difference. 
In accordance with this, we also find that the number of 
stars and cluster mass remaining at a particular time do not 
vary noticeably with the metallicity of the cluster. 

ISchulman et~ai1 l|2012l ) evolved TV-body models start- 
ing with TV = 8 192 stars and different metallicities to find a 
size difference between metal-poor and metal-rich clusters, 
in terms of the half-mass radius. This is in disagreement 
with our result s and those of the Monte Ca rlo models of 
iDownind (|2012t ). The ISchulman et al.1 l|2012t) models were 
evolved with some softening so that the effects of close bi- 
naries were not included. They were evolved to a dynamical 
age of 5£ r h which translated to physical ages in the range 
of 100 - 500 Myr for the small- TV models. The claim is that 
the results should be applicable to larger clusters, including 
GCs, because the impact of different stellar evolution and 
mass-loss histories at various Z will not depend on TV, and 
also because they performed models in the range of 1 024 
to 16 384 stars that showed similar half-mass radius evolu- 
tion. We would counter that as the MS lifetime of a MS 
turn-off star changes with age and the half-mass relaxation 
timescale of a cluster varies with TV, it is not at all obvious 
that the interplay between stellar evolution and cluster dy- 
namics will scale in a straightforward manner . Indeed, our 
model s here and the open cluster models of Irlurlev et al.l 
(2004) with TV ~ 30 000, both show that the half-mass ra- 
dius of metal-rich models can be smaller than that of the 
metal-poor models at early times (see Fig. [6]) but that the 
difference is erased or even reversed later in the evolution. 
Factors including different core-collapse times, the stellar 
evolution of low-mass stars as a function of metallicity (par- 



ticularly for globular clusters with ages of 10 Gyr or more) 
and different remnant masses need to be taken into account 
to gain the full picture. Furthermore, statistical fluctuations 
are generally prevalent in small-TV simulations and it can 
be necessary to average th e results of many in stances to es- 
tablish true behavior (e.g. iKiipper et al.|[200ct ). Our models 
presented here are at the lower edge of the GC mass function 
but even for these we would suggest that larger models again 
are desired before making any final judgment about the size 
measurements of GCs in general. Howeve r, our agreemen t 
with the large-scale Monte Carlo models of lDownind l|2012l ). 
performed with 5 x 10 5 stars, on the issue of half-mass radius 
variation (or non-variation) with metallicity is reassuring. 

In contrast to the evolution of the half-mass radius, we 
find that the half- light (or effective) radius does vary with 
metallicity. We find that blue, metal-poor clusters can ap- 
pear on average 17% larger than red, metal-rich clusters, 
with even larger differences possible when comparing indi- 
vidual models. This is in agreement with observations of 
extra-galactic GC systems, where size d i fferences of 17—30% 
llLarsen et al.ll2~00ll ; Ijordan et al. I l2005l ; IWoodlev fc Gomel 
l20ld ) have been found. It is also in agreement with the 
Monte Carlo models of lDowninel(|2012l ). Indeed, our TV-body 
models and these Monte Carlo models provide excellent in- 
dependent validation of the main result - that the observed 
size differences in GCs are likely caused by the interplay 
of stellar evolution and mass segregation. Stellar evolution 
causes low-Z stars to be brighter than their high-zT coun- 
terparts while mass segregation causes the most massive 
remnants to sink to the centre. Successively more massive 
remnants in \ow-Z clusters leads to a steeper surface bright- 
ness profile for high-zT clusters. The overall mass segrega- 
tion is similar for metal-poor and metal-rich clusters but 
more effective in the luminous stars for high-zT clusters ow- 
ing to a higher main-sequence turnoff mass. Th i s is in excel- 
lent agreement with the predictions of Ijorda n (2004) using 
multi-mass Mitchie-King models to estimate the size differ- 
ence between blue and red GCs, finding a difference of 14% 
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due to the combined effect of mass-segregation and stellar 
evolution. 

The apparent size difference does have a dependence 
on the treatment of remnants. When ejecting all NSs and 
BHs, no significant size difference (half-light radius) is found, 
partly owing to the fact that one of the variations with 
metallicity (remnant masses) has been negated. When we 
retain ~ 5% of the NSs and BHs arising from the primor- 
dial popjilaUon 2 _our results are in general agreement with 
the lDowninel (|2012h models that retained large numbers of 
BHs. While there are uncertainties in the retention frac- 
tions for NSs and BHs, there are also uncertainties for the 
masses of remnant BHs. We have us ed the stellar evolu- 
tion wing mass-loss prescriptions from iHurlev et all (2000), 
while improved, Z dep endent mass-loss rates are now avail- 
able l|Vink et al.ll200lh . However, the resulting differences 
for BH masses are mo st apparent for stars above 40 Mg 
(Belc zynski et af]|201Cf ). while just a few stars are drawn 
from this mass range in the models presented here. 

The average size difference of 17% implies that blue 
GCs do indeed appear larger as a result of metallicity effects. 
Since this is at the lower end of what is found in observa- 
tions, other causes (such as projection effects) can also be 
expected to play a role. In the future we plan to extend our 
study by performing additional iV-body simulations that ex- 
plore parameters such as larger N, smaller initial size and 
differing initial density profiles, as well as different cluster 
orbits, to further understand the effects of cluster evolution 
and environment on measured sizes. Our spread of individ- 
ual measurements in Fig.[8]can be compared to extragalactic 
studies of GC systems as well as in the Milky Way, in which 
half-light ra dii of GCs are found to be dis tributed between 1 
to 8 pc f e.g.lLarsen fe Brodidliool Fig. 4, ISpitler et afll2006l 
Fig. 19. iMadrid et al.l 120091 Fig. 10). Since clusters of dif- 
ferent masses and at different galactocentric distances are 
included in the observational samples, a larger scatter is 
expected than for our models (which currently give values 
between 2 — 6pc). We would expect the model spread to 
increase when we extend our study to include a range of 
cluster parameters. 

In addition to the half-light radius, we have also ana- 
lyzed the evolution of the mass-to-light ratio. When com- 
paring cluster models evolved purely through stellar (but 
no dynamical) evolution with the thorough iV-body mod- 
els, there is little change in M/L. As seen before in 
iBaumgardt fe Makinol (2003), we find that M/L increases 
with time, where dynamical interactions lead to a decrease 
in M/L as low-mass stars (carrying a high mass-to-light ra- 
tio) are preferentially lost from the cluster. The decrease in 
overall cluster luminosity with time results in an increase of 
the mass-to-light ratio. 
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